Statistical analysis of air and sea temperature anomalies 
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This paper presents a global air and sea temperature anomalies analysis based upon a combination 
of the wavelet multiresolution analysis and the scaling analysis methods of a time series. The 
wavelet multiresolution analysis decomposes the two temperature signals on a scale-by-scale basis. 
The scale-by-scale smooth and detail curves are compared and the correlation coefficients between 
each couple of correspondent sets of data evaluated. The scaling analysis is based upon the study 
of the spreading and the entropy of the diffusion generated by the temperature signals. Therefore, 
we jointly adopt two distinct methods: the Diffusion Entropy Analysis (DEA) and the Standard 
Deviation Analysis (SDA). The joint use of these two methods allows us to establish with more 
confidence the nature of the signals, as well as their scaling, and it yields the discovery of a slight 
Levy component in the two temperature data sets. Finally, the DEA and SDA are used to study 
the wavelet residuals of the two temperature anomalies. The temporal regions of persistence and 
antipersistence of the signals are determined and the non-stationary effect of the 10-11 year solar 
cycle upon the temperature is studied. The temperature monthly data cover the period from 1860 
to 2000 A.D.E. 
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I. INTRODUCTION 

The statistical analysis of time series is a challenging 
problem of statistical mechanics. This is due to the fact 
that there are still many unsettled problems. The most 
important seems to be that the techniques of analysis 
that are currently used are based on the assumption that 
the time series under study are generated by stationary 
processes. In general this is not the case. The time series 
mirroring complex processes are usually non-stationary 
in nature. The non-stationary condition seems to be a 
very general property, although it has any number of pos- 
sible sources in any system. For instance the origin of 
non-stationarity in the case of solar flares is given by the 
solar cycles (for a recent review about this interesting is- 
sue, see Ref. 0) and a special caution must be adopted 
to take the effects of this non-stationarity into account 
In fact, it has been recently shown || that the mem- 
ory left after detrending annual periodicity is much less 
intense than imagined in earlier publications [Q . Another 
issue, which seems to be still unsettled, is as to the sta- 
tistical nature of the fluctuations, once their genuinely 
stationary nature has been assessed. Are these fluctua- 
tions Gaussian? Are these fluctuation of Levy in nature? 

In this paper we want to illustrate an efficient approach 
to the solution of these difficulties. To stress the efficiency 
of this approach we apply it to the analysis of global air 
and sea temperature anomalies, a problem where, as we 
shall see, properly detrending non-stationary components 
is an essential request to shed light into the nature of the 
process under study. The approach we intend to use rests 



on the joint use of the Diffusion Entropy Analysis (DEA) 
and wavelet analysis of time series. DEA was born as 
an efficient way to detect scaling HHD, with applica- 
tions to sociological || and astrophysical (2) processes. 
This technique of analysis has been applied with suc- 
cess also to the study of DNA sequences [@,|| and heart 
beat rhythms in cardiac patients ||. Furthermore, some 
attention has been devoted to establish the connection 



between DEA and the Kolmogorov complexity |10| and 
it is becoming clear that this technique can also be used 
to study the transition from dynamics to thermodynam- 
ics, a crucial property that is used with success to study 
small portions of large sequences Q , thereby establishing 
a possible way to address the problem of non-stationality. 
Research work is currently being done to make it possible 
to utilize this technique to address the cases of multiple 
scaling Jill . 

Wavelet techniques are a powerful method of analysis 
]l2[ that localizes a signal simultaneously in time and fre- 
quency. We use wavelets for the purpose to decompose 
the signal in smooth, detail and residual components. 
The wavelet decomposition has been shown to be an ef- 
ficient way of detrending from the data a non-stationary 
component in a natural way, so as to bypass the main dif- 
ficulties concerning the non-stationary nature of the data 
under study jl3| . The adoption of DEA makes the scal- 
ing emerge and also sheds light into the statistical nature 
of the fluctuations around the non-stationary bias. 

Let us now illustrate the time series under study in this 
paper. The time series of annually averaged global sur- 
face temperature anomalies have attracted the attention 
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of many scientists since the pioneering work of Nicolis 
and Nicolis juj . Dynamical systems theory has provided 
a new quantitative perspective on the predictability of 
weather and climate processes. For more recent attempts 
along these lines the interested reader can consult the 



work of Rcf. |15|. The conceptual paradigm behind our 
analysis is that of intermittence as a dynamic source of 
Levy statistics jy] . The detection of Levy scaling Q af- 
ter detrending the non-stationary component would cor- 
roborate the validity of an intermittent dynamic model. 

It is worth remarking that the term temperature 
anomalies is a technical definition adopted in the current 
literature on weather and climate processes on earth to 
denote air and sea temperature departures from a mean 
temperature value. Therefore, this term must not be con- 
fused with the term anomalous diffusion that is related 
to our conceptual paradigm. We use the data on tem- 
perature fluctuations to generate a diffusion process that 
is compared to the standard Brownian motion. The de- 
parture of this resulting diffusion process from Brownian 
diffusion is called anomalous diffusion. The DEA aims at 
measuring the strength of this anomaly. Consequently, 
we can say that one of the aims of this paper is to deter- 
mine the anomalous nature of the diffusion process gen- 
erated by air and sea temperature anomalies. The data 
analyzed are updated continuously by the Climate Re- 
search Unit in the United Kingdom [|l6| and the Hadley 
Centre for Climate Prediction and Research, Meteoro- 
logical Office |l7j ]. The earliest of attempts at collecting 
these data were done in 1986 |1| and has culminated 
in what is recognized as one of the most accurate data 
files for global air temperature, and global sea surface 
temperature (SST) p9| . The land-air temperature has 
been corrected for non-climatic errors, such as changes 
in the location of the weather stations and changes in 
instrumentation E(J . The SST data have been corrected 
for changes in instrumentation that was used before 1942 
pl|-^3| . The data consists of numbers that represent the 
departure from a mean temperature in order to see a 
change from a global average of sorts. The mean tem- 
perature used was the 1961-1990 mean temperature. 

Figs. 1 show the global air (a) and sea (b) temperature 
anomalies in the period from 1860 to 2000 A.D.E. The 
dashed lines are the wavelet multiresolution S7 smooth 
portion of the signal. The exact mathematical definition 
of wavelet smooth curve is given in Section 3. To un- 
derstand the meaning of the data illustrated in Fig. 1, 
it is enough for the reader to consider the S7 smooth 
as a convenient way to establish a type of mean value 
about which the temperature fluctuations take place. In 
fact, as it will become clear in Section 3, the S7 smooth 
is obtained by a wavelet average of the data over time 
intervals of 128 months. 

The global air and sea temperature anomalies are very 
similar to one another. The mean value of temperature 
fluctuation shows that from 1860 to 1915 the average 



temperature is almost constant, with a change of only 
0.2° Celsius. From 1915 to 1950 the temperature in- 
creases by 0.6° Celsius. From 1945 to 1980 the aver- 
age temperature remains almost constant again. Finally, 
from 1980 to current time there is a further increase of 
the average temperature of almost 0.4° Celsius. What 
about the fluctuations around this mean value? Figs. 2 
show the spectral density against the period in months 
of the global air (a) and sea (b) temperature anomalies 
in the period from 1860 to 2000 A.D.E. . Some of these 
periodicities are reported in Table I. The main periodic- 
ities involve a time period 6 and 12 months long (related 
to the yearly cycle of the Earth orbiting the Sun), 9-12 
years long, 21-22 years long and, finally, a strong pe- 
riodicity 55-57 years long (the last three cycles are all 
established solar cycles with the 9-12 year cycle being 
the most widely known and observed one). 

This paper aims at settling several questions concern- 
ing these data. The correlation among air and sea tem- 
perature smooth curves seems to be evident. However, it 
is not so clear if these correlations exist also at the level 
of fluctuations. Is there anomalous scaling? If there is 
anomalous scaling, does it rest on smooth curves or fluc- 
tuations? Which is the statistical nature of these fluctu- 
ations? Are they Gaussian fluctuations? Are they Levy 
fluctuations? 

The outline of the paper is as follows. In Section 2, to 
make the paper as self contained as possible, we make a 
short review of wavelet multiresolution analysis. Section 
3 is devoted to the study of the correlation between global 
air and sea temperature anomalies at the scale of both 
the wavelet smooth and detail curves. Section 4 is a short 
review of the DEA and shows this technique at work on 
the data without any decomposition. Section 5 illustrates 
the joint use of wavelets and DEA. Finally, in Section 6 
we make a balance on the results obtained in this paper. 



II. WAVELET MULTIRESOLUTION ANALYSIS 

Wavelet analysis (l2) is a powerful method to analyze 
time series that is attracting the attention of an ever 
increasing number of investigators. Wavelet Transform 
makes use of scaling functions, the wavelets, which are 
characterized by the important property of being local- 
ized in both time and frequency. These functions inte- 
grate to zero and, usually, are normalized. A scaling 
coefficient t characterizes a wavelet. The length 2t mea- 
sures the width of the wavelet and defines the time scale 
analyzed by the wavelet. Two typical wavelet functions 
that are widely used in the continuous wavelet transform 
are the Haar wavelet and the Mexican hat wavelet Jl2| . 
The Haar wavelet is defined as: 
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-1/a/2t, t - t < u < t 
i> Tt t(u)={ l/\/2r, t<u<t + r 
0, otherwise 



(1) 



The Mexican hat wavelet is the second derivative of a 
Gaussian function. Given a signal the Continuous 

Wavelet Transform is defined by 



W(j,t) = / iJ T , t (u) £(u) du 



(2) 



The original signal can be recovered from its Continuous 
Wavelet Transform via 
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W(r,t) i> T ,t(u) dt 



where is a constant that depends on the wavalet func- 
tion ]l2| ]. The double integral of Eq. (j|) suggests that 
the original signal may be decomposed in "continuous 
details" that depend on the scale coefficient r. How- 
ever, there exists a discrete version of the wavelet trans- 
form, the Maximum Overlap Discrete Wavelet Transform 
(MODWT), which is the basic tool needed for studying 
time series of N data via wavelet. In the Ref. p2| , the 
reader can find all of the mathematical details. For the 
purpose of this paper, it is important to have in mind 
only one of the important properties of the MODWT: 
the Wavelet Multiresolution Analysis (WMA). It is pos- 
sible to prove that given an integer Jo such that 2 J ° < N, 
where TV is the number of the data points, the original 
time series represented by the vector X can be decom- 
posed as follows: 



Jo 



with 



Si -i 



(4) 



(5) 



The detail Dj represents changes on a scale of 2r = 2 J , 
while the smooth Sj represents averages on a scale of 
2tj = 2 Jo . We term wavelet residuals the quantities 
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Jo 



X-5> 



Jo 



(6) 



It is then evident that we can interpret the residuals as 
fluctuations about the local mean value evaluated on the 
time scale 2tj — 2 Jo . At this stage, the reader should 
fully understand the comments made in Section 1 about 
the data illustrated in Fig. 1. 

For the reader to properly appreciate the value of 
WMA, we show this technique at work by means of the 



results illustrated in Figs. 3 and 4. These figures show 
the WMA of the global air and sea temperature anoma- 
lies in the years 1860-2000. The analysis is done by us- 
ing the Daubechies least asymmetric scaling wavelet fil- 
ter (LA8) Q. The LA8 wavelets look similar to the 
Mexican hat but they are asymmetric, a fact that makes 
them more plastic than the Mexican hat wavelet. We 
have plotted the WMA for Jo = 7. Figs. 3 compare the 
smooth curves S4, S5, S6 and S7 and the detail curves 
D4, D5, D6, D7 of the air (solid lines) and sea (dashed 
lines) temperature data. Figs. 4 show the details Dl, 
D2, D3 and D4 of the two sets of data. Details Dj show 
the fluctuations of the temperature on a scale of 2r = 2 J 
months. According to Eqs. (f|) and (j^), the sum of all 
seven details and the 57 smooth give the two original 
signals. The 56 smooth curve is given by 57 + D7; the 
55 smooth is given by 56 + D6, and so on, till to 51, 
given by 52 + D2. 

Figs. 3 show that the global air and sea temperature 
anomalies are closely correlated to each other but they do 
not coincide. If it happens that in a given time region air 
is hotter than sea, immediately afterward the opposite ef- 
fect takes place, and sea is warmer than air. The smooth 
curves of Figs. 3, and especially the S7 and S6 smooth 
curves, show that in the time regions of persistent tem- 
perature increase (time regions 1860-1880, 1910-1950 and 
1980-2000) air is hotter. Instead, in the time regions of 
persistent temperature decrease (time regions 1880-1910 
and 1950-1980) air is colder. This may be explained by 
the fact that heat capacity of water is higher than heat ca- 
pacity of air, thereby implying that it takes more time for 
water temperature either to increase or to decrease. The 
D4, D5, D6, D7 details illustrate the temperature fluc- 
tuations of both air and sea corresponding to the time 
scales of 16, 32, 64 and 128 months respectively. The 
fluctuations of both kinds of data, air and sea, look re- 
markably similar. The analysis of the details Dl, D2 and 
D3 of Figs. 4 show, instead, that the fluctuations of the 
air temperature are larger than those of the sea temper- 
ature. This effect, too, may be related to the higher heat 
capacity of water. Finally, Figs. 4 show a stronger fluc- 
tuation of the data during the period 1860-1880. There 
is also a strong fluctuation of sea temperature in the D3 
detail, in the 1920-1950 time region. Further study is re- 
quired to assess whether these stronger fluctuations are 
due to natural phenomena or to some artifact of a non 
conventional way of data acquisition. 



III. MULTIRESOLUTION CORRELATION 
ANALYSIS 

The use of the Multiresolution Correlation Analysis via 
wavelet is a simple procedure [Q. We decompose the 
two temperature datasets into rests and details by using 
WMA as demonstrated in the previous section. Then, 
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we create pairs of partners, the first component of the 
pair being either a rest or a detail of the sea temperature 
data and the second the corresponding rest or detail of 
the air temperature data. For any given pair of datasets 
(xi,yi): i = l,...,N, the linear correlation coefficient r is 
given by the formula 

V // E l fa -x? \/Ei {yi~y? 

where, as usual, x is the mean of the former sequence 
and, y is the mean of the latter sequence. The value 
of r lies between -1 and 1, and can include the extreme 
values 1 and —1. It takes on a value of 1, termed "com- 
pletely positive correlation," when the data points lie on 
a perfect straight line with positive slope, with x and y 
increasing together. The value 1 holds independently of 
the magnitude of the slope. If the data points lie on a 
perfect straight line with negative slope, y decreasing as 
x increases, then r has the value -1; this is called "com- 
pletely negative correlation." A value of r near zero in- 
dicates that the variables x and y are uncorrelated. The 
two sequences are considered to be significantly corre- 
lated when the value of \r\ is close to 1. 

Fig. 5 illustrates the correlation coefficient, r, between 
air and sea temperature anomalies as a function of the 
wavelet index j, and thus, for the reasons illustrated in 
Section 3, as a function of the wavelet scale 2r = 2- 7 . The 
top curve denotes the correlation between the air and the 
sea Sj smooth curves, with j ranging from to 10. The 
bottom curve denotes the coefficient of correlation be- 
tween air and sea details, Dj, with j ranging from 1 to 
10. The value r = 0.87 corresponding to the top curve 
at j = is the correlation coefficient between the two 
original temperature data without any filtering. For the 
sake of reader's convenience the values of r are reported 
in Table II. From Fig. 5 we see the small details Dl, 
D2, D3 and D4 are not significantly correlated. These 
details refer to a time scale until 16 months, and conse- 
quently we can conclude that within this time range there 
are no significant correlations between sea and temper- 
ature fluctuations. The correlation between the details 
increases with increasing the scale index and becomes sig- 
nificant for D5, D6 and D7. The reader can consult Figs. 
3, which confirms in fact visually that these details are 
significantly correlated. The correlation between the D7 
details yields a local maximum, r — 0.95. All this means 
that the S4-S7 smooth curves are the best indicators of 
the correlations between air and sea temperature anoma- 
lies. In fact, according to the prescriptions of Section 3, 
the S4-S7 smooth curves are obtained by adding to S4 the 
D5, D6 and D7 details, respectively, and these details, 
as shown by the bottom curve of Fig. 5 are correlated. 
Again, the reader can make a visual inspection of these 
smooth curves by consulting Figs 3. Details smaller than 
the D5 detail have small correlations. This also means 



that to shed light on the properties that make sea tem- 
peratures different from air temperatures we must focus 
on small details. Smooth curves larger than S7 are very 
well correlated but do not afford additional information, 
because they are too smooth. 

IV. SCALING ANALYSIS 

Scale invariance has been found to hold empirically 
for a number of complex systems p5[ ] and the correct 
evaluation of the scaling exponents is of fundamental im- 
portance to assess if universality classes exist |2^| . A 
widely used method of analysis of complexity rests on 
the assessment of the scaling exponent of the diffusion 
process generated by a time series. See, for instance, 
Refs. H||f7],|]|27) . According to the prescription of Ref. 
|p7f , we interpret the numbers of a time series as generat- 
ing diffusion fluctuations and we shift our attention from 
the time series to the probability distribution function 
(pdf) p(x,t), where x denotes the variable collecting the 
fluctuations. The scaling property takes on the form 

f(*.*) = ^(J) ' ( g ) 
where 6 is the scaling exponent. 

A. Gauss and Levy diffusion 

There are two main forms of anomalous diffusion. The 
first is the generalization of Brownian motion, proposed 
years ago by Mandelbrot [^5| , known as Fractional Brow- 
nian Motion (FBM) and yielding for the diffusion process 
a variance increasing in time as t 2H . This kind of anoma- 
lous diffusion fits the scaling definition of Eq.(^) with 
8 = H and F(y) being a Gaussian function of y. A second 
form of anomalous diffusion is obtained by generalizing 
the Central Limit Theorem (CLT). The prescriptions of 
the Generalized Central Limit Theorem (GCLT) J2f| are 
as follows. Let us assume that the diffusing variable x 
is the sum of t independent random variables £i, each of 
which has a probability distribution, symmetric around 
£ = 0. Let us assume also that for large values of |£|, this 
distribution is an inverse power law, with index ji > 1, 
so as to fit the normalization condition, and fi < 2, so 
as to violate the CLT constraint. Then, according to the 
GCLT, for t — ► oo , the diffusion becomes stable and the 
Fourier transform of p{x, t) gets the form 

p(k,t)=exp(-b\k\ a t), (9) 

where & is a kind of generalized diffusion coefficient, de- 
termined by the strength of the fluctuations and a — 
/i + 1 . It has been shown [^8| that there is scaling. This 
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means that the resulting diffusion process fits the condi- 
tion of Eq. (0), with 5 = 1/a, namely, 



1 



(10) 



Both forms of anomalous diffusion are idealization of 
reality. The theory of FBM implies that the correspond- 
ing scaling property holds true at any time scale, while its 
dynamic derivation |29[| , in the persistent case H > 0.5, 
suggests that it is a time asymptotic property generated 
by the fluctuations of a Gaussian variable whose 
correlation function has inverse power law with index 
(3 < 1 (H = 1-/3/2). The main problem with this 
dynamic interpretation is that there are no stable theo- 
rems behind the Gaussian property of the "microscopic" 
fluctuation £(£). It is necessary to supplement the dy- 
namic model with random ingredients that do not have 
a dynamic origin pjJ. The earlier illustrated approach 
to Levy statistics, called Levy flight, is judged to be un- 
realistic, because it involves a random velocity that can 
be arbitrarily large. To bypass this difficulty recourse 
is given to Levy walk pl| . We adopt the following dy- 
namic model jnj. Let us consider a sequence {ri,Sj}, 
with i — 0,1,..., co. The numbers are random num- 
bers with the distribution density 



V>(r) = (/i - 1) 



{T + tY 



(11) 



where T is a positive constant. Note that to ensure the 
stationary condition the additional condition /i > 2 is 
necessary, and we make this assumption also in this pa- 
per. Thus, we have /i < 3 to ensure the anomalous char- 
acter of the resulting diffusion process and [i > 2 to make 
our dynamic picture stationary. It is easy to prove that 
< t >= T/(fi — 2). The numbers Si have the values 
1 and —1, determined by the coin tossing rule. For a 
generic time t, let us consider the time tjq fitting the con- 
ditions ijv = T o + T i + ■■■ t n-i +t n < t and tjy + t n > t. 
Then the trajectory prescribed by this dynamic model is 
x(t) = W[r s + Tisi + ...T N ^ 1 s N -i + (t-t N )s N ]. We see 
that in this case, due to the large memory time, the ran- 
dom walkers can travel ahead or backwards by quantities 
with the same distribution as the Levy flight. However, 
this takes a time proportional to the traveled length. The 
correlation functions of each elementary jumps, either W 
or — W has an inverse power law with index /3, which is 
now given by (3 = \x — 2. It is important to notice that 
the process is now multiscaling [jllj , due to the fact that 
the propagation fronts propagate linearly in time. The 
central part of the distribution, if we neglect the trun- 
cation produced by the finite velocity of the propagation 
front, is given by 8 = l/(/x — 1), namely, by Eq.(|10|). 

It is important to stress that this dynamic model to 
Levy statistics yields, at any finite time t, finite second 
moments. Consequently, it could be interpreted as a form 



of FBM. However, in this case H would not correspond to 
the correct scaling of the central part of the distribution 
and 5 and \x are related to this pseudo-scaling H by M . 



8 = 



1 



and 



respectively. 



3- 2H 



2H, 



(12) 



(13) 



B. The diffusion algorithm 



Let us consider a sequence of N numbers 



1. 



(14) 



The goal is to establish the possible existence of a scaling, 
either normal or anomalous, in the most efficient way 
as possible without altering the data with any form of 
detrending. First of all, let us select an integer number 
t, fitting the condition 1 < t < N. This integer number 
will be referred to by us as "diffusion time" . For any 
given time t we can find M(t) = N — t + 1 sub-sequences 
defined by 



Si — Si+sj 



ith s = 0, 



,N — t. 



(15) 



For any of these sub-sequences we build up a diffusion 
trajectory, s, defined by the position 



(16) 



1=1 



The direct evaluation of variance is probably the most 
natural method of variance detection. All trajectories 
start from the origin x(t = 0) = 0. With increasing time 
t, the sub-sequences generate a diffusion process. At each 
time t, it is possible to calculate the Standard Deviation 
of the position of the M(t) sub-sequences with the well 
known expression: 



D(t) = 



M(t) - 1 



(17) 



where x(t) is the average of the positions of the M(t) 
sub-trajectories at time t. The exponent H is defined by 



D(t) oc t H 



(18) 



We call this approach to the scaling evaluation Stan- 
dard Deviation Analysis (SDA) . In Ref . |j| the interested 
reader can find an illustration of the traditional tech- 
niques of scaling detection and of why all of them are 
virtually equivalent to the SDA. 
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The DEA, is based upon the following algorithm. We 
have to partition the cc-axis into cells of size e(t). When 
this partition is made, we have to label the cells. We 
count how many particles are found in the same cell at 
a given time t. We denote this number by Ni(t). Then 
we use this number to determine the probability that a 
particle can be found in the i-th cell at time t, Pi(t), by 
means of 



(19) 



At this stage the entropy of the diffusion process at the 
time t is determined and reads 



£(*) = -!>(*) HMt)} ■ 



(20) 



The easiest way to proceed with the choice of the cell size, 
e(t), is to assume it to be a fraction of the square root 
of the variance of the fluctuations £(i), and consequently 
independent of t. If the scaling condition of Eq. (|J) holds 
true, it is easy to prove that 



S(t) =A + 8 In(t) 



(21) 



where, in the continuous approximation, 

/oo 
dy F(y) ln[F(y)] , (22) 
-oo 

with y = x/t 5 . The scaling Eqs. ( [i"8| ) and ( |2l] ) determine 
the exponents H and S. 



N = 1680 data, a number with a square root of about 40. 
The pictures show that for t = 50 the standard deviation 
is almost 40 times larger than the standard deviation at 
the first step of diffusion. It seems, therefore, that the 
statistics are rich enough to get a satisfactory pdf and 
consequently reliable scaling properties. 

The high values of the exponents imply a strong per- 
sistence. This means that the temperature changes grad- 
ually month by month. The fact that the exponents 
stemming from the sea temperature anomalies are higher 
than those produced by the air temperature anomalies 
means that the sea temperature anomalies are character- 
ized by a persistence higher than that of the air tempera- 
ture anomalies. This can be explained as an effect of the 
higher heat capacity of the water. These results confirm 
the results of Sec. 3, shown in Figs. 4, where the fluc- 
tuations of the air temperature anomalies at short scale 
are stronger than those of the sea temperature anoma- 
lies. Finally, we note that the both exponents H a , for air, 
and H s , for sea, are larger than S a and S s respectively. 
This means that the pdf of the two diffusion processes 
is a little bit larger than a Gaussian distribution. The 
four exponents fulfill the Levy Walk Diffusion relation 
dl2| ) within the accuracy of our statistical analysis. This 
means that global air and sea temperature anomalies not 
only are characterized by Levy statistics, but are a mani- 
festation of the dynamic model illustrated in Section 4 A. 
This means that the alternated periods Tj of high and low 
temperature are distributed according an inverse power 
law with fi < 3. According to Eqs. ( p"3] ) and (|To|) , we ob- 
tain fi a = 2.13 ± 0.02 and fj, s = 2.08 ± 0.02. The results 
are summarized in Table III. 



C. Data analysis 

Fig. 6 shows the numerical results by using the DEA 
(Figs. 6a and 6c) and the SDA (Figs. 6b and 6d) of 
the global air (Figs. 6a and 6b) and sea (Figs. 6c and 
6d) temperature anomalies. In the ordinate axis we plot 
D(t)/D(l) and 5(t) - 5(1). Thus, the curves start from 
1 and 0, respectively. The straight lines are function of 
the type foE(t) = 8 m W an d fso(t) = t H and become 
straight lines as a consequence of the linear-log (DEA 
case) and log-log (SDA case) representations we are 
adopting. The global air temperature anomalies are char- 
acterized by a pdf scaling coefficient 8 a = 0.87 ±0.02 and 
a standard deviation scaling coefficient H a — 0.92 ±0.01. 
The global sea temperature anomalies are characterized 
by a pdf scaling coefficient S s — 0.89 ± 0.02 and a stan- 
dard deviation scaling coefficient H s = 0.94 ± 0.01. The 
figures are plotted for a period t = 50 months. For value 
of t larger than 50, saturation effects due to the statisti- 
cal property appear. The fitting is a kind of mean result 
obtained by averaging the results corresponding to fitting 
the first 20, 30 and 40 points. The upper time limit con- 
dition of t — 50 is dictated by the fact that we are using 



V. WAVELET MULTIRESOLUTION DIFFUSION 
ANALYSIS 

In this section we introduce a method of analysis, based 
on the joint use of wavelet decomposition and the diffu- 
sion approach to scaling, the latter method resting on 
both DEA and SDA. This method turns out to be pow- 
erful, and we refer to it as Wavelet Multiresolution Dif- 
fusion Analysis, WMDA. Figs. 7 and 8 show the SDA 
and the DEA of the global air (a) and sea (b) tempera- 
ture applied to the residuals Rj, where, as in the earlier 
sections, j indicates the wavelet scale index. Each resid- 
ual contains all details at smaller scales. Therefore, the 
WMDA allows us to determine the diffusion spreading 
at each time scale, as stemming from the corresponding 
details, an important piece of information. The residuals 
Rj are obtained by detrending the original data with the 
smooth curves Sj obtained with the wavelet multireso- 
lution analysis; see Eq. (^). In each of the four figures 
there are nine curves. The figures (a) refer to air, the 
figures (b) refer to sea. Figs. 7 refer to SDA and Figs. 8 
to DEA. From top to bottom the curves of these figures 
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denote the original data (1), and the residuals (2) R9 (2) 
, (3) R8 (3) (9) R2(9) . The curves of Figs. 7 and 
8 are an illustration of WMDA and prove that this is a 
useful tool to study complexity of a dynamical system. 
Let us analyze them in detail. 

a- Figs. 7 and 8 look similar but are not identical. The 
curves afforded by DEA are more detailed than those 
obtained by using SDA. This is because DEA is more 
sensitive to fluctuations than standard deviation, which 
only measures the spreading of the diffusion trajectories. 

b- The straight lines correspond to the function 
fsD(t) = t - 5 , in Figs. 7, and to the function foE(t) = 
0.5 In(t), in Figs. 8, and look straight due to the adoption 
of the log-log (Figs. 7) and linear-log (Figs. 8) represen- 
tations. These straight lines serve the purpose of signal- 
ing what would be the behavior of an ordinary Gaus- 
sian diffusion. According to the Mandelbrot interpreta- 
tion , the curve with slopes larger than that of these 
straight lines indicate persistent diffusion or supcrdiffu- 
sion, and the curves with slopes smaller than that of these 
straight lines indicate antipersistent diffusion or a subdif- 
fusion region. By comparing the curves of Figs. 7, resting 
on SDA, to the corresponding curves of Figs. 8, resting 
on DEA, we see that the times at which the SDA curves 
cross the line fsoif) — t°' 5 are larger than the times at 
which the DEA curves cross the line fDEif) = 0.51n(t). 
This is a further proof that the dynamics behind the data 
cannot be Fractional Brownian Motion. This kind of dif- 
fusion process would imply SDA and DEA crossing the 
lines denoting ordinary diffusion at the same time. The 
figures suggest that the diffusion is characterized by a pdf 
whose tails are more persistent than the Gaussian tails 
of Fractional Brownian Motion. The analysis of Section 
5C proved that the diffusion is anomalous and fulfills the 
Levy Walk condition. In this section, with the help of 
WMDA we show that the data are characterized by an 
anomalous dynamics at short as well as long time scales. 

c- By comparing the curves of Figs. 7a and 8a to the 
corresponding curves of Figs. 7b and 8b, we note the air 
temperature curves always have a slope smaller, even if 
slightly smaller, than the slope of the corresponding sea 
temperature curves. This means that at each wavelet 
scale the higher heat capacity of the water makes the sea 
temperature data more persistent than the air ones. This 
is related to the higher heat capacity of the water. 

d- Figs. 7 and 8 show that each curve is characterized 
by one leading periodicity. The dynamical reasons for 
this property are easily accounted for by noticing that 
a given periodicity of the data causes a periodic conver- 
gence of distinct trajectories. After an initial spreading, 
with a consequent increase of both variance and entropy, 
there are incomplete regressions of the initial condition. 
Since each curve corresponds to a given scale, the ob- 
served processes of regression correspond to the leading 
periodicity of that temporal scale. Any residuals Rj con- 
tains all details at smaller scales, and, as a consequence, 



the leading periodicity is not necessarily related to the 
wavelet temporal scale Tj by simple relations. For ex- 
ample, the R7 and R8 curves of Figs. 7a and 8a show 
almost the same periodicity but T7 is one half of rg. More- 
over, the main periodicities do not coincide exactly for 
the global air and sea temperature anomalies. The fig- 
ures show that the year periodicity and its multiples have 
a strong effect until the residual R6. A periodicity of 17- 
20 years dominates the R7 and R8 residuals. Finally, a 
periodicity of 57-59 years characterizes the R9 rest. Figs. 
7 and 8 show that WMDA may be an interesting comple- 
ment to the spectral density analysis of Figs. 2 because 
it shows the main periodicity and its contribution to the 
information that characterizes that scale level. 

e- Figs. 7 and 8 show that both entropy and standard 
deviation of any residual converges to a horizontal line. 
This is due to the detrending of the smooth part of the 
data, Sj, which makes the hidden periodicities show up. 
At the same time, these hidden periodicities imply the 
deterministic nature of the signal and consequently yield 
entropy saturation. In other words, the fluctuations & of 
the residuals data Rj can generate trajectories with only 
a limited spreading. The height of the horizontal lines 
measures the maximum spreading (in the case of SDA, 
Figs. 7) and the information or entropy (in the case 
of DEA, Figs. 8) for each time scale. It is interesting 
to notice that the shorter the time scale the faster the 
transition to saturation. This means that the role of 
periodicities becomes more and more important as we 
decrease the time scale. 

Table IV summarizes some of the information con- 
tained in Figs. 7 and 8, and can help the reader to un- 
derstand the balance of the results of this paper that will 
be given in Sec. 6. 

Finally, the reader may wonder what happens if we de- 
trend one of the wavelet details from the original data. In 
our analysis we observe that the scaling properties coin- 
cide, within the limit of the analysis statistical accuracy, 
with those of Section 4 C for t < 50. This applies to 
detrending from the original data any detail between Dl 
and D10. This means that the wavelet smooth curves 
S determine the scaling because they determine the per- 
sistent properties of the signal. The scaling for t < 50 
is not conditioned by any cycle in the data. If we con- 
sider interval larger than t = 50 we observe the largest 
discrepancy by detrending the detail D7. Figs 9a (air) 
and 9b (sea) show the results. The detrended data have 
a lower entropy for 30 < t < 200. This entropy seems to 
scale with 8 close to 0.5 until t = 600, like a Brownian 
diffusion. However, the lack of sufficiently rich statistics 
in that region does not allow us to get any conclusion 
about the real value of the scaling. The SDA, instead, 
does not detect this difference, Figs 9c and 9d; the en- 
tropic analysis turns out to be more sensitive than the 
variance analysis. The increasing of entropy of the orig- 
inal data for 30 < t < 200 have a non-stationary origin 
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due to a cycle that is detrended by the detail D7. In 
fact, detail D7 corresponds to the wavelet time scale of 
128 months, that is, 10-11 years, and it contains the im- 
portant 10-11 years solar cycle. 



VI. CONCLUSION 

The key results of these paper can be summarized as 
follows: 

(i) Understanding the various aspects of a global cli- 
mate system. Global climate is a complex system in 
which many factors affect one another and without a de- 
tailed look into the behavior of the Earth different tem- 
perature systems and their correlations, any attempt to 
understand it will be difficult at best. For example, un- 
derstanding these behaviors is important to comprehend 
local and global climate changes, as well as, to attempt 
to rebuild past climate more accurately |33j |. In fact, our 
analysis shows that the differences between the behavior 
of the data on the shorter time scale lends itself to cau- 
tion when performing these reconstructions. The under- 
standing and reconstruction of climate to times previous 
to instrumental records will allow us to gain some knowl- 
edge of the long term behavior of future climate |34|. We 
show, in fact, that no significant correlation exists for de- 
tail up to time intervals of 16 months. For larger time 
intervals this correlation becomes visible. Our analysis 
has shown the differences in the fundamental types of 
statistical behavior of both regions of the Earth. This 
knowledge could be used to further our ability to recon- 
struct past climate in an attempt to better understand 
our dynamic global environment. 

(ii) Scaling and its relevance for the air and sea temper- 
ature anomalies. The joint use of DEA and SDA allowed 
us to establish better the nature of the signals, as well as 
their scaling. The data seems to have a slight Levy com- 
ponent. The higher scaling of sea signal is interpreted 
as a consequence of the higher thermal capacity of sea 
water. 

(iii) Periodicity effects on scaling. We proved that scal- 
ing detected by using the DEA is not affected by periodic- 
ities. The smooth curves are responsible for scaling and 
scaling is not influenced by cycles. The smooth curves 
are responsible for a steady increase of entropy of the 
diffusion. 

(iv) Residuals and periodicities. We showed that peri- 
odicities emerge at the level of residuals, namely the por- 
tions of the signals obtained detrending the smooth parts. 
The diffusion entropy of residuals saturates, thereby im- 
plying that after a given time there is no further infor- 
mation increase. We can therefore, conclude that there 
should be no concern about a possible influence of period- 
icities on scaling. The scaling found is a genuine property 
and we can freely adopt it as an indicator of correlations. 



The sea temperature data, yielding a higher scaling, im- 
ply a larger correlation and the statistical analysis of this 
paper makes compelling this important conclusion. 

(v) Details and non- stationary effect. We show that 
detail D7, that contains the important 10-11 year solar 
cycle, causes a non-stationary effect that is detected by 
DEA but not by SDA. This shows that DEA has a higher 
sensitivity. 

(vi) Joint use of entropy and decomposition. We think 
the benefit of the joint use of DEA and wavelets is evi- 
dent. The wavelet decomposition generates a set of new 
time series, corresponding to tuning the wavelet micro- 
scope to a given time scale, and the DEA establishes the 
information of these components, and makes it evident 
why periodicities set an upper limit on entropy increase. 

It is interesting to point out that the earlier analysis of 
temperature anomalies @,|l5| has been done using deter- 
ministic chaos and the evaluation of the Lyapunov coef- 
ficient, and so the Kolmogorov-Sinai (KS) entropy. Our 
research work seems to support the dynamical model of 
Ref. |ll| , which is a form of Levy walk. This means a 
connection with turbulence and intermittence. The val- 
ues of y, emerging from this analysis are very close to 
(i = 2 and thus to the border with the non-stationary 
dominion |35| ] . It is known that the KS entropy vanishes 
for (j, < 2. The authors of Ref. Q seem to rest on a 
condition of vanishing KS entropy to address the intrigu- 
ing problem of climate predictability. We think therefore 
that the joint use of DEA and Compression algorithms 
JToj l , applied after the wise detrending method illustrated 
in this paper, might contribute further processes towards 
the ambitious goal of predictability. 
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Figure 1 
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FIG. 1. Global air and sea temperature anomalies in Celsius degree (years: 1860-2000). Fhe dashed lines are the wavelet 
multiresolution smooth curves S7. 
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Figure 2 
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FIG. 2. Global air (a) and sea (b) temperature anomalies spectral density analysis against the period. In Table I there are 
some of the main periods. 
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TABLE I. Main periods present in the global air and sea temperature anomalies. The values are in years. 
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Figure 3 
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FIG. 3. Wavelet Multiresolution Analysis of the global air (solid lines) and sea (dashed lines) temperature anomalies in 
Celsius degree. The figures show the smooth curves S7, S6, S5, S4 and the details curves D7, D6, D5 and D4 that are 
associated to the scales of 256, 128, 64 and 32 months. 
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Figure 4 
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FIG. 4. Wavelet Multiresolution Analysis of the global air (left) and sea (right) temperature anomalies in Celsius degree. 
The figures show the details D4, D3, D2, Dl that are associated to the scales of 16, 8, 4, 2 months. 
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Figure 5 
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FIG. 5. Correlation coefficient r between the global air and sea temperature anomalies against the wavelet scale index. The 
top curve denotes the correlation coefficient between the wavelet smooth curves, from SO to S10. The smooth curves SO refer to 
the original data without any filtering. The bottom curve denotes the correlation coefficient between the wavelet detail curves 
from Df to DIO. See Table II for the value of r. 
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TABLE II. Correlation coefficient r between the global air and sea temperature anomalies against the wavelet scale index. 



The value r — 0.87 of SO is the correlation coefficient between the global air and sea temperature anomalies without any 
filtering. 
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Figure 6 
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FIG. 6. DEA (Figs, a and c) and SDA (Figs, b and d) of the global air (Figs, a and b) and sea (Figs, c and d) temperature 
anomalies. The straight lines correspond to functions of the type ]de{£) = 5 ln(t) and fso(t) = t H , which become straight 
lines in the linear-log representation of this figure. The global air temperature anomalies are characterized by a pdf scaling 
coefficient 8 a = 0.87±0.02 and a standard deviation scaling coefficient H a = 0.92±0.01. The global sea temperature anomalies 
are characterized by a pdf scaling coefficient S s = 0.89 ± 0.02 and a standard deviation scaling coefficient H a = 0.94 ± 0.01. 
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TABLE III. Exponents H, 8 and fi for the global air and sea temperature anomalies. 
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Figure 8 
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FIG. 8. DEA of the global air (a) and sea (b) temperature anomalies of the residuals From top to down, the curves 
are the SDA of (1) the original data, (2) R9, (3) R8,..., (9) R2. The straight line is fDE(t) = 0.51n(i) that corresponds to the 
Gaussian diffusion. 
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TABLE IV. Summary of the information contained in Figs. 7 and 8. The first four columns report the values of the standard 
deviation and of the entropy at the first step of diffusion for both global air and sea temperature anomalies and their rests Rj 
for j = 2, 9. The four following columns, from the 5th to the 8th, report the values of the height of the horizontal lines that 
measures the maximum spreading (in the case of SDA, Figs. 7) and the information or entropy (in the case of DEA, Figs. 8) 
that corresponds to each wavelet scale. The reported heights are relative to the values of the SDA and DEA at I = 1. This 
means that the SDA and DEA heights are defined by Dn ne = Dj/D(l) and Sum — Sj — S(l), respectively, where Sj and Dj 
are the values of the entropy and of the standard deviation of the diffusion process generated by the residuals Rj . The last two 
columns report the main periodicities present in each residuals Rj . 
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Figure 9 
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FIG. 9. DEA and SDA of the global air (a, c) and sea (b, d) temperature anomalies of the original data (solid line) and 
of the data detrended of the detail D7 (dashed line). The straight lines in (a) e (b) have a slope of 5 = 0.51 ± 0.02 (a) and 
S — 0.53 ± 0.02 (b) that correspond to the Gaussian diffusion. The effect of the detail D7 is detected by DEA but not by SDA. 
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